Phase and frequency entrainment in locally coupled phase oscillators with repulsive interactions 
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Recent experiments in one and two-dimensional microfluidic arrays of droplets containing Belousov - 
Zhabotinsky reactants show a rich variety of spatial patterns [J. Phys. Chem. Lett. 1, 1241-1246 (2010)]. 
The dominant coupling between these droplets is inhibitory. Motivated by this experimental system, we study 
repulsively coupled Kuramoto oscillators with nearest neighbor interactions, on a linear chain as well as a ring 
in one dimension, and on a triangular lattice in two dimensions. In one dimension, we show using linear sta- 
bility analysis as well as numerical study, that the stable phase patterns depend on the geometry of the lattice. 
We show that a transition to the ordered state does not exist in the thermodynamic limit. In two dimensions, 
we show that the geometry of the lattice constrains the phase difference between two neighboring oscillators to 
2n/3. We report the existence of domains with either clockwise or anti-clockwise helicity, leading to defects 
in the lattice. We study the time dependence of these domains and show that at large coupling strengths the 
domains freeze due to frequency synchronization. Signatures of the above phenomena can be seen in the spatial 
correlation functions. 

PACS numbers: 05.45.Xt,89.75.-k,82.40.Bj 



I. INTRODUCTION 

Synchronization JJ, 0] , in which individual oscillators col- 
lectively organize into some phase and frequency relationship 
with each other is prevalent in nature. It has been studied ex- 
tensively in many biological and chemical systems. Examples 
of systems which exhibit partial or complete synchronization 
include oscillations in neuronal networks yH7|]> and coupled 
chemical oscillators 181- Hlll . 

The simplest model that exhibits synchronization was pro- 
posed by Kuramoto [12], and applies to weakly-coupled sys- 
tems where the phase varies slowly compared to the intrin- 
sic frequencies. This is an exactly solvable model. Even 
though the model is mean-field (every oscillator coupled to 
every other oscillator), and ignores amplitude variations lfl3ll . 
it has been very successful in capturing the qualitative features 
of synchronization. Locally coupled versions of the Kuramoto 
model have been studied using numerical techniques, and by 
mapping to statistical mechanics models Il3l4l8ll . 

Most of the well-studied models invoke attractive (or ex- 
citatory) coupling between the oscillators. In these systems, 
a transition from ordered to disordered phases is seen at fi- 
nite coupling strengths, and the lower critical dimension for 
frequency entrainment was shown to be d = 2 liMl El- 
More recently, oscillators with repulsive (or inhibitory) cou- 
plings ha ve g enerated interest due to applications to biological 
systems i20H23ll . Unlike the attractive coupling case, where 
all oscillators relax into an in phase state irrespective of the ge- 
ometry, the solutions in the case of oscillators with inhibitory 
coupling depend strongly on the underlying geometry, since 
the connectivity of the lattice can frustrate the local order pre- 
ferred by the coupling. In a one dimensional chain, the na- 
ture of the global pattern was seen to change when the cou- 
pling was changed from a local to a global coupling l24ll . The 
global inhibitory coupling is "frustrated" since every oscilla- 
tor prefers to be n out of phase with every other oscillator. 
Other variants of repulsively coupled systems have also been 
studied l25ll . Synchronous, traveling waves which undergo a 



transition to spatiotemporal chaos were seen in coupled circle 
maps with repulsive coupling f26ll . The role of frustration in 
determining patterns and tuning adaptive networks has been 
studied in the context of biological systems |27H30I1 . 

Recently, an experimental system with inhibitory coupling, 
and controllable geometry and coupling strength, was real- 
ized in an array of water microdroplets surrounded by an 
oil medium IBll433ll . These water droplets contained reac- 
tants of the oscillatory Belousov-Zhabotinsky (BZ) reaction. 
Bromine, which is a constituent reactant, dissolves preferably 
in the oil medium, and provides the inhibitory coupling be- 
tween the droplets. The strength of the coupling between the 
droplets is varied by varying the size of the droplets. In one di- 
mensional array of these droplets, anti-phase synchronization 
as well as Turing patterns are observed. In two dimension ge- 
ometry, these droplets rearrange into a hexagonal array. At 
weak coupling strengths, the droplets relax into a "27r/3" pat- 
tern in which any two neighboring droplets are 120° out of 
phase. At large coupling strengths, a "7T-S" pattern is seen, 
in which the central droplet relaxes into a non-oscillatory sta- 
tionary state, and the surrounding six droplets in the hexago- 
nal array exhibit anti-phase oscillations. This experiment was 
qualitatively explained by including diffusion of inhibitory 
species in the FKN model generally used to describe BZ reac- 
tion in. 

Motivated by this experimental system, we study a local 
variant of the Kuramoto model with repulsive coupling. We 
show that in one dimension (ID), the neighboring oscillators 
prefer to relax into an anti-phase state. The spatial pattern of 
phase differences change when the geometry is changed from 
a linear chain to a ring. For the ID ring, the phase pattern de- 
pends on the number of oscillators. We explain the observed 
patterns by studying attractors at infinite coupling using linear 
stability analysis. We show that there is no synchronization 
transition in the thermodynamic limit since the critical cou- 
pling strength scales with the number of oscillators. In two 
dimensions (2D), we study oscillators on a triangular lattice, 
where we expect the effects of frustration to be the most pro- 
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nounced. We show that the oscillators prefer to relax to the 
27r/3 state seen in the experimental system. We report the 
existence of domains with either clockwise or anti-clockwise 
helicity, in which phases of any three neighboring oscilla- 
tors either increase or decrease in a given direction, leading 
to domain-wall defects in the lattice. We study the time de- 
pendence of these domains and show that at large coupling 
strengths the domains freeze into a "glassy state" in which the 
phases are disordered but the frequencies are synchronized. 
We characterize these phenomena by studying spatial correla- 
tion functions. Finally, we discuss these results in the context 
of the experimental BZ microdroplet system. 

The paper is organized as follows. We give details of the 
model in section|II] We discuss the results from our ID study 
in section [III] In section IIV1 we discuss the results obtained 
in 2D. We conclude with a discussion of our results in the 
context of the experimental system. 




FIG. 1 : (Color online) Order parameter, A*, for N = 3(0), N = 7(0). 
and N = 10(A) for linear chain averaged over random initial phase 
configurations. 



II. THE MODEL 



III. RESULTS IN ID 



The BZ micro-oscillators were modeled using locally cou- 
pled Kuramoto oscillators placed on a lattice. The equations 
governing the phase of the oscillators are: 



-ojj + K^ sin(0, - <pj) 

(ij) 



(1) 



Here, the oscillator at site i is coupled locally to its nearest 
neighbors {j}. The intrinsic frequency of the individual os- 
cillators is given by a>j. This frequency is chosen randomly 
from a Gaussian distribution with zero mean and variance 
cr, and is a source of quenched disorder in this system lfl3ll . 
The strength of the coupling between the sites is given by K. 
As discussed below, the effective coupling strength that con- 
trols the synchronization behavior is K/cr, and therefore, we 
set cr to unity without loss of generality. When the coupling 
strength in the above system is attractive, viz. K < 0, the 
oscillators synchronize in phase with each other. However, 
when K > 0, the coupling is repulsive and resulting attractors 
depend intrinsically on the underlying geometry of the lattice. 

In ID, we study a linear chain, as well as a ring geometry. 
In 2D, we study a triangular lattice, with six nearest neighbor 
connections. Both in ID and 2D, the linear size of the array 
is L = 64 unless otherwise specified, with the lattice spacing 
taken to be unity. The numerical integrations are performed 
using Runge-Kutta algorithm with a fixed step-size 0.05. We 
have checked the step-size dependence of our results, and ex- 
cept in the glassy state, which will be discussed in detail be- 
low, the results are not sensitive to reducing the step size be- 
yond this value. Statistical properties of the data are obtained 
by averaging over 100 realizations of quenched frequencies in 
ID and over 50 realizations in 2D. 



A. Linear Stability 

In ID we consider the system with two different bound- 
ary conditions: free ends (linear chain), and periodic (ring). 
With the free end boundary conditions, the oscillators at ei- 
ther end of the array have only one neighbor to which they 
are coupled, whereas the oscillators form a ring in the case of 
periodic boundary conditions, in which the first oscillator is 
coupled to the N' h . As we will show, this can create frustra- 
tion in the system, and give rise to interesting spatial patterns 
when the oscillators are repulsively coupled. 

We begin our study in ID by performing a linear stability 
analysis of the deterministic version of Eq. Q] In this system, 
there is no disorder in the intrinsic frequencies: cr = 0, and all 
the oscillators have the same frequency, u>. Since the effective 
coupling strength controlling phase and frequency patterns is 
K/cr, taking the limit of cr — > is, therefore, equivalent to 
taking \K\ — » oo. The linear stability analysis is thus a study 
of the stability of the attractors in the infinite coupling strength 
limit. 

The linear stability analysis is most transparent if we per- 
form a change of variable to the phase difference = 0, - <p j+l 
between nearest neighbors. The transformed equations in ID 
are: 



9/ = -A'sin(0 / _i) + 2Ksm(Gi) -Ksm(6 i+l ). 



(2) 



It is known 1 34] that in a linear chain, Eqnf2]has 2 N ~ l attrac- 
tors, but only one stable solution which depends on the sign of 
K. For K < only the completely in-phase solution, 6>, = 0, 
is stable, while for K > only the completely anti-phase so- 
lution, 9j = n is stable. 

It is under periodic boundary conditions that we find a much 
richer long time behavior. Solving for the fixed points, we find 
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the condition 

sin(0*) = sin(6»*) (3) 

for all i, / Additionally, we have the constraint imposed by 
the boundaries 

N 

J] & = 2m (4) 

;=1 

where n is an integer. Thus the system has fixed point solu- 
tions for all values 8* = 2nn /N. We can determine the stabil- 
ity of these solutions by noting that the Jacobian is a circulant 
matrix OSll . for which, the eigenvalues are easily computable. 
Specifically, we find 

A,„ = 4Ksin 2 (mn/N)cos8*; m = 0, 1, ...,N- 1. (5) 

Since sm 2 (rrm/n) is always positive, the stability depends only 
on the phase difference 8* , and the sign of K. For values of 
K > 0, all values of 8* between n/2 and 3n/2 are stable, where 
as for K < 0, 8* between -n/2 and n/2 are stable. It should be 
noted that while 6* — n is the most stable phase configuration 
for K > 0, it is not accessible for systems with odd numbers 
of oscillators. Contrast this with the case of K < where 
the most stable state, 8* = 0, is accessible for any number of 
oscillators. 

We studied the synchronization transition in ID by numer- 
ical integration of the dynamical equations for the phase (Eq 
|2J. Here, we take the intrinsic frequencies of our oscillators 
to be chosen randomly from a Gaussian distribution. The mean 
of the Gaussian is often set to zero, but the behavior of the sys- 
tem is the same for any value of the mean. There are two dis- 
tinct processes leading to the complete synchronization of the 
system: frequency entrainment, where the oscillators adjust 
their frequency of oscillation until all evolve with the same 
frequency, and phase ordering, where the oscillators develop 
a uniform phase difference between all neighboring elements. 
One interesting question that we address is whether the fre- 
quency and phase synchronization occur as two separate tran- 
sitions, or is there a single transition where both phases and 
frequencies lock in. 
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FIG. 2: (Color online) Order parameter, A2n/3, for N = 3 on a ring. 
Three different initial phase configurations are shown: all initially 
in-phase(O), all initially —2n/3 out of phase(A), and initially n/2 out 
of phase(<>). 
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FIG. 3: Bifurcation diagram of time averaged oscillator frequencies 
<f>i plotted as a function of coupling strength K for a ring of 64 oscil- 
lators. 



B. Phase Ordering 

To investigate the phase ordering we define a phase order 
parameter fl4il . as 

1 N 

A 9 = <-| jy^-U-^D (6) 
J'=l 

which is unity for the fully ordered system, and zero when 
the system is completely disordered. Here, q is the expected 
phase difference between oscillators and the average is over 
realizations of w,. 

As discussed above, we expect the system with free end 
boundary conditions to always go to the completely anti-phase 



state at high enough coupling strength. In Fig. Q]the order pa- 
rameter, is plotted as a function of coupling strength for 
different system sizes. Although the initial phase configura- 
tions were chosen at random, A n still tends to a value of one 
with increasing K, showing that there is only one attractor, 
8* = n. Additionally, we can see that as is increased a 
larger coupling is required to reach the completely ordered 
state. This agrees with the conclusions reached in Ref. Il4ll 
regarding attractively coupled oscillators, that that there can 
be no complete phase synchronization in the thermodynamic 
limit for d < 5 111 . 

With periodic boundary conditions, the story becomes a bit 
different. Fig. [2]shows the order parameter, A2 K /j,, for a three 
oscillator ring. The linear stability analysis of this system 
gives two stable states, 8* = +2n/3. We notice in the fig- 
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FIG. 4: (Color online) Critical coupling strength as a function of 1/N 
where N is the number of oscillators, for linear chain(O) and ring(O) 
geometries. 



ure that the behavior of the order parameter depends upon the 
initial phases of the oscillators. If all oscillators are given the 
same initial phase, 0q — where 6q = </>, - fa-i, they are 
equally likely to fall into either of the +2n/3 attractors, thus 
takes on a value of roughly one half. If on the other 
hand, < 0q < n the system always goes to the 2n/3 state and 
A2^/3 tends to one, while for n < 6q < 2n, the system goes 
to the -2n/3 state and A2n/3 tends to zero. These two states 
in the three-oscillator ring appear in the 2D, triangular lattice 
as the two helicities of the 2n/3 state, and lead to the frozen 
domains, as discussed in detail in section HVl 



C. Frequency Entrainment 

When the coupling between the oscillators is turned off, 
each of the oscillators in the system evolves in time according 
to its own frequency, w,. As we switch on and increase the 
coupling, we observe the oscillators to form local frequency 
entrained clusters. As the coupling is increased further, the 
clusters merge with one another until eventually, at some crit- 
ical value K c , all of the oscillators are entrained with a com- 
mon frequency. This clustering process is illustrated in Fig. [3] 
By carrying out our study at different values of the random- 
ness, <x, we have directly tested that the parameter controlling 
the synchronization is K/cr: increasing the variance of the w, 
distribution has the same effect as reducing the coupling. As 
Fig. |4]shows, the critical coupling strength obtained from fre- 
quency synchronization increases as the number of oscillators 
is increased, and there is no finite-/^ transition in the thermo- 
dynamic limit. We note that linear chain requires a larger cou- 
pling strength to obtain frequency entrainment than a ring sys- 
tem with the same number of oscillators and K c shows no sys- 
tematic dependence on odd/even values of N for either geom- 
etry. Clustering behavior in the frequency synchronization of 
a ring of oscillators has been observed earlier in a locally cou- 
pled Kuramoto model with attractive coupling |36ll . In gen- 




FIG. 5: Phase difference plots(left) and space-time plots(right) at 
four different coupling strengths for the ring geometry. The space- 
time plots show a selection of 15 oscillators taken from a 64 oscilla- 
tor system. Time is along the horizontal axes and space is along the 
vertical axes. 



eral, we find that even though the character of frequency syn- 
chronization in this repulsively coupled system is qualitatively 
similar to that seen in attractively coupled systems, the phase 
patterns obtained can be significantly different. 

Fig. [3] does not give any information about how the os- 
cillators in a given cluster are distributed in space. Fig. [5] 
provides insight into how the clustering of frequencies occurs 
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FIG. 6: (a),(c) Time-series of the phase tf>(t) of a typical oscillator in the interior and at the boundary of a domain respectively, at coupling 
strength K = 3.0. (b),(d) Power spectrum |/ v | 2 calculated from the above time-series. The power spectrum \f v \ 2 shows a power-law decay with 
an exponent ~ 1.7. The inset to the plots (b),(d) show the corresponding auto-correlation function A(t) The solutions were sampled after an 
initial transient t = 1000 over a time interval T = 10000. 



spatially. The figure shows four space-time plots over a range 
of coupling and a corresponding plot of the phase difference 
as a function of position on the lattice. At low coupling both 
the frequencies and phases are disordered. There is, however, 
a structure visible in the space-time plots even at K/cr — 1 
showing that the oscillators in the same frequency cluster are 
spatially correlated. As the coupling is increased to K/cr = 10 
we can see that all the oscillators evolve with approximately 
the same frequency, but maintain some disorder in the phase 
relationships. It is not until a coupling strength two orders 
of magnitude larger that we see perfect phase ordering at this 
system size. The patterns that we obtain from our model is 
qualitatively similar to that observed in experiments on BZ 
droplets JmIH. 



at different locations in the lattice, at a coupling strength K = 
3.0, is shown (As will be discussed later, these oscillators are 
located in the interior and at the boundary of a phase-locked 
domain respectively). The two oscillators exhibit widely dif- 
ferent dynamics. While the first oscillator shows relatively 
monotonic behavior, the latter oscillator exhibits intermittent 
dynamics in which the oscillator follows a monotonic trajec- 
tory in between bursts of chaotic dynamics. The power spec- 
trum |/ v | 2 of these oscillators obtained by taking a Fourier 
transform of the time series </>(t) (f v — J dt cp(t) exp(/vf)), ex- 
hibits a power-law decay |/ v | 2 ~ with an exponent ( ~ 1 .7 
(Fig. |6{b),(d)). This spectrum indicates the existence of a 
multitude of time-scales in the phase dynamics of the repul- 
sively coupled oscillators. The exponent ( varies between 



IV. RESULTS IN TWO DIMENSIONS 

In the two dimensional experimental setup lO , the 
droplets containing the BZ reactants reorganize into an hexag- 
onal array. We model this situation by placing the Kuramoto 
oscillators on a triangular lattice, each oscillator being repul- 
sively coupled to its six nearest neighbors. As was shown in 
the one-dimensional model, the phases of the oscillators tend 
to align n out of phase with their neighbors, when they are 
repulsively coupled. The locally-preferred order of n phase 
difference between nearest-neighbor oscillators is not glob- 
ally realizable on the triangular lattice which has fundamen- 
tal loops containing three oscillators that are mutual nearest 
neighbors; a situation analogous to the three-oscillator ring 
discussed in section [III] As we show below, this frustration 
leads to (a) much richer dynamics in comparison to attrac- 
tive coupling, (b) defects in the phase order, and (c) a phase- 
disordered yet frequency synchronized state. 

The rich and spatially heterogenous dynamics in the case 
of repulsively coupled oscillators is evident in Fig. |6ja),(c), 
in which the time-series of the phases <p(t) of two oscillators 
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FTG. 7: Distribution of the magnitude of phase differences between 
two adjoining oscillators on a row P(A<f>), plotted for various cou- 
pling strengths K. The distribution is a gaussian which peaks around 
a phase difference of 2n/3. The inset shows the rescaled distribution 
P(Atp)K- a plotted against (A0 - 2n/3)K a . a = 2/3. 
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FIG. 8: (Color online) (a) Graphic illustrating the different helicities seen on the triangular lattice, (b) Domains of opposite helicities repre- 
sented by the red (darker) and white regions seen for coupling strength K = 3.0. The green (lighter) regions indicate the frequency entrained 
oscillators, (c) Colorscale map of the phase differences between two oscillators on a row. The dark and lighter regions represent location of 
oscillators with a phase difference of -2n/3 and 2n/3 respectively. 



1.7 - 2.0 for different oscillators, and is typical of the spatial 
heterogeneity seen in the lattice. The auto-correlation func- 
tion A(r) calculated by taking an inverse Fourier transform of 
the power-spectrum |/ v | 2 , A(t) = J dv\f v \ 2 exp(-i'vT), shown 
in the insets to Fig. |6jb),(d), shows a slow relaxation. 

Non-exponential relaxations and multiple time scales are 
a hallmark of glassy systems 13711 . and reflect the presence of 
many metastable states, and a complex free-energy landscape. 
In the coupled oscillators system, the multiple time scales pre- 
sumably arise from a complex attractor landscape. Within the 
phase-coupled model, we therefore expect to see glassy be- 
havior, and the correlation functions discussed below, indeed 
indicate frozen, disordered phase patterns. 

At weak coupling strengths K, the phases of the oscillators 
are mainly disordered, though a few neighboring oscillators 
show a tendency to oscillate n out of phase with each other. 
As the coupling strength is increased above K > K p 1 .5, the 
oscillators relax into a state in which each oscillator is 2n/3 
out of phase with its neighbor. This resembles the 2n/3 state 
seen in the BZ micro-oscillators setup 113 311 . The distribution 
of the magnitude of the phase differences between neighbor- 
ing oscillators P{A<p), shows a peak at 2tt/3 (Fig. |7}. The 
distributions at different coupling strengths can be collapsed 
onto a single curve, when the phase differences are rescaled 
as given below. 

P(A4>, K) = K- a f((Acf> - 2n/3)K a ) (7) 

The inset to Fig. [7] in which the rescaled distribution of phase 
differences is plotted for different coupling strengths, illus- 
trates the scaling behavior. The form of the scaling function 
f{{A<p - 2n/3)K") is seen to be gaussian. Hence, the variance 
of the distributions is proportional to l/K 2 ", where a = 2/3. 
This suggests that as the coupling strength approaches infin- 
ity, the distribution approaches a (5-function centered at 27r/3. 
Hence, oscillators prefer to relax into a state with a phase dif- 
ference of 27r/3. 



Interestingly, however, at coupling strengths K > K p , do- 
mains with opposite helicities form in the lattice. In each 
domain, the phases of any three neighboring oscillators vary 
continuously in either clockwise or an anti-clockwise direc- 
tion. The two opposite helicities have been illustrated by a 
graphic in Fig. [Ha). The two different domains on the lattice, 
have been represented by the red (darker) and white regions 
in Fig. [8jb). The phase differences between two neighbors in 
a row in the two domains could be either 27r/3 or -2jt/3 de- 
pending on the helicity. This is shown in Fig. |Uc), where the 
phase differences in the lattice has been plotted. The darker 
regions represent location of oscillators with a phase differ- 
ence of -2n/3 between the neighbors. At the boundaries of 
these domains, the phases of the neighboring oscillators are 
seen to be aligned anti-parallel with each other. So, in a lat- 
tice, a non-zero fraction of neighboring oscillators are n out 
of phase with each other, and the fraction of n out-of-phase 
neighbors is proportional to the fraction of sites belonging to 
the domain boundaries. In the interior of a domain, the helic- 
ity of the phases are frozen. However, the oscillators at the do- 
main boundaries can fluctuate between either of the helicities, 
since they exist in an anti-phase state. This is evident in the 
time-series of the phase <p(t) of a boundary oscillator plotted 
in Fig |6jc), in which the phase exhibits intermittent dynam- 
ics when the domain boundary shifts positions in the lattice. 
These domain boundaries then act as defects in the lattice. At 
weak coupling strengths (K p < K < 3.5), the domains show a 
tendency to coarsen with time. In other words, domains with 
one of the helicities grow resulting in a shrinkage in the size of 
domains with the opposite helicity. The selection of the dom- 
inant helicity is unbiased in a set of random initial conditions. 
As the coupling strength is increased further (K > 3.5), the 
domains freeze after a very short transient, and do not change 
with time. 

In analogy with Ising spin systems i38h . we define a helic- 
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FIG. 9: Helicity order parameter m plotted as a function of the cou- 
pling strength K. The error bars represent the fluctuations in the order 
parameter (Am 2 ). The data has been collected after a transient time 
t = 5000. 



ity order parameter that characterizes the domains on the lat- 
tice. We measure the phase relation between any three neigh- 
boring sites 1,2,3 as shown in the graphic in Fig. [8{a), and 
check the helicity of this triplet. We associate a helicity index 
p(R) = +1 or -1 with site 1 if the helicity is clockwise or anti- 
clockwise respectively. If N + and N- are the total number of 
sites in each domain, we define the helicity order parameter 
m as m — \N+ - N-\/N. The dependence of the helicity order 
parameter m on the coupling strength is shown in Fig. [9] The 
helicity order parameter m is close to zero when the phases 
are disordered. As the coupling strength K is increased, do- 
mains with opposite helicities are formed. Due to the coars- 
ening of the domains and the subsequent domination of one 
of the helicities on the lattice, the helicity order parameter m 
approaches 1. On further increasing the coupling strength, 
frozen domains with finite sizes are seen, and an order param- 
eter m < 1 is obtained. A clearer signature of the onset of the 
domain formation is seen in the fluctuation of the helicity or- 
der parameter. When the phases are disordered, viz K < K p , 
the fluctuation in the order parameter {Am 2 ) = {{m - m) 2 ) is 
negligible. When distinct domains start forming in the lattice, 
as the coupling strength is increased beyond K p , the fluctua- 
tion in the order parameter {Am 2 } shows a distinct increase, as 
shown in Fig. [9] 

Signatures of the dynamics of domains can be detected in 
the time-dependent spatial correlation function C(r, t), which 
is defined as C(r, t) = J dRp(R, t)p{R+r, t), where p(R, t) is the 
helicity index +1,-1 associated with the site at position/?. At 
weak coupling strengths (K p < K < 3.5), where the domains 
coarsen with time, the correlation function C(r, f) shows an 
increasing correlation length as time progresses (Fig. [TUta)). 
The domain size L(t) scales with time t as L(t) ~ f, where 
j] = 0.5. This is seen in the inset to Fig. [TOf a). in which the 
correlation function obeys the scaling form C(r, t) ~ f(r/P). 
The scaling function f(x) is seen to be an exponential of the 
form f(x) = exp(-ajt). An exponent of r\ — 0.5 is usually 
seen in the case of non-conserved scalar fields, where the do- 
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FIG. 10: (a) Spatial correlation function C{r,t) plot- 
ted at coupling strength K = 3.0 and at times t = 
1024(+),2048(x), 4096(A), 8192(a), 12288(V), 163840,20480(0). 
The inset shows the correlation function C(r, t) plotted against the 
scaled variable r/f, where 77 = 0.5. The dashed line represents 
an exponential fit to the data, exp(-3.4x). (b) Spatial correlation 
function C(r,t) plotted for coupling strength K = 5.0 at times 
t = 4096(A), 8192(x), 12288(V), 16384(a), 20480(0). 
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main growth is surface tension-driven l39ll . At larger coupling 
strengths, the correlation function does not evolve with time 
after a transient period, concurrent with the frozen domains 
observed at these parameter values (Fig. fTUl b)). 

The structure factor S (k, f), defined as the fourier transform 
of the spatial correlation function C(r, f) corroborates the for- 
mation of distinct domains in the lattice. As shown in Fig. 
[TT1 the structure factor S(k,t) decays as S{k,t) ~ k~ e ,9 
2.74 + 0.04 ^ 3 for large wavevectors k. This is in accor- 
dance with Porod's law wherein the structure factor decays 
as S(k,t) ~\Jk d+l whenever there are well-defined domain 
boundaries [40]. 

A possible reason for the coarsening of domains at weak 
coupling strengths can be inferred from the collective behav- 
ior of the frequencies of oscillators. We define an effective 
frequency, to e fj for each oscillator as 

U eff = (4>i(to + T) - 4>i(t ))/T (8) 

Here, the frequencies have been averaged over the interval t = 
5000 after discarding data for a transient time to = 1000. The 
inset to Fig. [T2l shows a bifurcation diagram, in which all the 
frequencies of the oscillators have been plotted as a function 
of the coupling strength K. At weak coupling strengths, the 
frequencies of the oscillators are not entrained, and a broad 
distribution of the oscillator frequencies is obtained. In this 
interval of coupling strength, the phases of the oscillators also 
evolve with different frequencies, hence we see the fluctuating 
domain boundaries, and growing domains. In this regime, the 
frequency entrained oscillators are seen in the interior of the 
domains, as shown in Fig. Ob). 

At larger coupling strengths, the oscillator frequencies are 
well entrained, and hence the phase patterns on the lattice are 
frozen. A quantitative measure of frequency entrainment is 
given by the frequency order parameter F a UM . which is de- 
fined as 

F a =< NJN > (9) 



where N„ is the maximum number of oscillators with identi- 
cal frequencies, and N is the total number of oscillators. When 
all the oscillators are frequency entrained, the order parame- 
ter approaches 1 . The frequency order parameter F w has been 
plotted as a function of coupling strength K in Fig. Q~2] We 
see that the coupling strength K = 3.5 at which complete fre- 
quency entrainment is seen, or F u = 1, coincides with the oc- 
currence of frozen domains. Hence, frequency entrainment is 
responsible for the frozen domains in this lattice at larger val- 
ues of the coupling strengths. The freezing of the domains is 
consistent with the form of the power spectra of phases shown 
in Fig. [6] which indicates a broad distribution of time scales 
in the system. The relative phases between neighboring oscil- 
lators can, therefore, be frozen in a disordered state while the 
whole system oscillates with a common frequency. The fun- 
damental reason for the appearance of a broad time-spectrum 
is not clear from our studies. We, however, expect that this 
feature is related to the tendency of neighboring oscillators to 
prefer a phase difference of n, thus freezing in domain bound- 
aries where this configuration is realizable. We are continuing 
to investigate the slow dynamics in this model. 



V. DISCUSSION 

In this paper, we studied locally coupled Kuramoto oscil- 
lators with repulsive coupling in one and two dimensions. 
In comparison to the attractively coupled system, the sys- 
tem studied shows much richer dynamics, with hints of mul- 
tiple time-scales, and a complex attractor landscape, which 
strongly depends on the underlying geometry of the lattice. 

In one dimension we showed that while the linear chain 
has one stable phase configuration, when periodic boundary 
conditions are introduced a range of attractors become avail- 
able to the system which depend upon the number of oscil- 
lators. The spatial patterns obtained qualitatively match with 
anti-phase patterns seen in the BZ micro-droplets experimen- 
tal setup. Additionally, we showed that even though the fre- 
quencies were entrained, the phase patterns were still disor- 
dered. Hence, the onset of frequency entrainment occurs at a 
lower coupling than the phase ordering. In two dimensions, 
we showed the existence of phase patterns similar to the 27r/3 
state seen in the BZ micro-oscillator system. We showed the 
existence of domains with clockwise and anti-clockwise he- 
licities in the same lattice. These domains showed coars- 
ening behavior at weak coupling strengths, where a grow- 
ing length scale could be detected that reached system size 
at large times. As the coupling strength was increased, we 
found that these domains freeze, such that the phase pattern 
does not change in time. This was attributed to frequency en- 
trainment at large coupling strengths, which ensured that the 
frozen phase patterns on the lattice oscillate with a common 
frequency. Mapping discrete dynamical systems to statisti- 
cal mechanics models gives new insights into the behavior of 
these systems fiH l42tl . The repulsively coupled Kuramoto 
oscillators serve as a good paradigm in which techniques and 
ideas related to statistical mechanics can be applied to an in- 
herently dynamical system. 
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Since Kuramoto phase oscillator formalism excludes the 
discussion of amplitude variations, we do not hope to see the 
7T-S state seen in the experimental system in this model. How- 
ever, we are looking at an extension of this model, which al- 
lows the oscillators to "switch" on or off based on its phase 
environment. 
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